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1. Introduction 

The PDWHEG method, first suggested in ref. [1], has been successfully applied to Z pair 
production [2], heavy-flavour production [3] and e + e~ annihilation into hadrons [4]. In 
ref. [5] a general description of the method was given, and in particular its implementation 
within the Catani-Seymour (CS) subtraction scheme [6] and within the Frixione-Kunszt- 
Signer (FKS) [7, 8] approach. 

In this paper we present an implementation of the W and Z hadroproduction cross 
section in the POWHEG framework, using the CS subtraction formalism. All next-to-leading- 
order (NLO) calculations used in POWHEG until now have been performed in the FKS 
method. In view of the popularity of the CS scheme, we find desirable to explore more in 
detail its use within POWHEG. In ref. [5] an outline of the implementation of the Drell-Yan 
production cross section in POWHEG in the CS scheme was given. In the present work we 
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depart slightly from that approach. In particular, we use a more appropriate form of the 
hardness variable used for the generation of radiation. As a further point, for the case of 
W production, if angular correlations in decay products are correctly taken into account, 
a new problem arises. In fact, the Born-level W cross section vanishes when the fermion 
decay products are exactly in the opposite direction of the incoming quark-antiquark pair, 
which causes a problem in the generation of radiation within the POWHEG method. We show 
that this problem has a simple solution, that can be easily generalized to all cases in which 
the Born cross section vanishes. 

The paper is organized as follows. In sec. 2 we describe how we performed the calcula- 
tion for the NLO W and Z cross section. In sec. 3 we discuss the POWHEG implementation 
and how to deal with vanishing Born cross sections. In sec. 4 we show our results for several 
kinematic variables and compare them with MCONLO [9] and PYTHIA 6.4 [10]. Finally, in 
sec. 5, we give our conclusions. 



2. Description of the calculation 

2.1 Kinematics 

2.1.1 Born kinematics 

We begin by considering the Born process for the annihilation of a quark and an antiquark 
into a lepton-antilepton pair 1 q + q — ► l + l. Following ref. [5] , we denote by k B and k e the 
incoming quark momenta, and by k\ and k 2 the outgoing fermion momenta. We call K 9 
and K e the incoming hadron momenta and define the momentum fractions x@ as 

k@ = x @ K @ . (2.1) 

We choose our reference frame with the z axis along the k e direction. We introduce the 
following variables 

a xM2 v 1 i (fci + fc 2 )° + (fci + k 2 f 

M = {h + h2) ' Y = 2 log ( kl+k2 r-(k 1+ k 2 r (2 - 2) 

that characterize the invariant mass and rapidity of the virtual vector boson. 2 We also 
introduce the angle 9i that represents the angle between the outgoing lepton and the k 9 
momentum, in the centre-of-mass frame of the lepton pair. The azimuthal orientation of 
the decay products is irrelevant here, since the cross sections do not depend upon it. We 
thus fix it to zero. At the end of the generation of the event, we perform a uniform, random 
azimuthal rotation of the whole event, in order to cover all final-state phase space. The 
set of variables M 2 , Y and 9i fully parametrize our Born kinematics. From them we can 
reconstruct 

X(B = \~s~ e ' Xe = \l~s~ e ' ( ' 



1 In case of W production the quark-antiquark and lepton-antilepton pairs have different flavour. We 
focus here for simplicity on leptonic decays of the vector bosons. Hadronic decays are treated similarly. 

2 The virtuality of the lepton pair M 2 will be distributed according to a Breit-Wigner formula around 
the squared mass of the vector boson My (where V stands for either the W ± or the Z). 
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where S = (K B + K e ) 2 . The leptons' momenta are first reconstructed in the longitudinal 
rest frame of the lepton pair, where each lepton has energy equal to M/2 and where the 
lepton momentum forms an angle Q\ with the e direction and has zero azimuth (i.e. it lies 
in the z,x plane and has positive x component). The leptons' momenta are then boosted 
with boost angle Y. 

The Born phase space in terms of these variables is written as 

(2.4) 

2.1.2 Real-emission kinematics 

The real emission process is described by the final-state momenta h\ , k 2 and k 3 , where k\ 
and &2 have the same meaning as before, and k 3 is the momentum of the radiated light 
parton. In the POWHEG framework, applied in the context of the CS subtraction method, 
one introduces a different real phase-space parametrization for each CS dipole. In the 
present case, we have two CS dipoles, with the two incoming partons playing the role of 
the emitter and the spectator. We consider the case of the © collinear direction. Thus, the 
emitter is the incoming parton with momentum k 9 . We introduce the variable 

g = 1 _ (fce + fce)-fr (25) 



k m ■ kr 



and the momenta 



K = ki + k 2 = k 9 + k e - k 3 (2.6) 
K = xk S) + k e . (2.7) 

Observe that K 2 = K 2 , which is the condition that fixes the value of x. When k% is 
collinear to fc we have 

xk 9 = fe e - k 3 , (2.8) 
and K = K. Following ref. [6], we introduce the boost tensor 

A W ^- «*ffi'^ + g££. (») 
the barred momenta 

k^ = K^ u {K,K)k u r r = l,2, (2.10) 

the barred-momentum fractions 

X ffi = IIj, Xq = x e , (2-11) 

and the barred incoming momenta 

fee = xk 9 = Xe-Kei k e = k e = x e K e . (2-12) 
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The barred momenta characterize the underlying-Born kinematics. We define then 

^ = (ll + = <* + w >, % - \ bg |±g^|±|£ , ( , 13) 

and the angle 0j is defined as in the Born case, but in term of the momenta k @ , k e , k\ and 

The radiation variables are given by 

v = ^£, <f>, (2-14) 

where is the azimuth of k$ around the z direction. 

From the set of variables M 2 , Y" e , x, v and we can reconstruct the full production 
kinematics for the real-emission cross section. We summarize the reconstruction proce- 
dure from ref. [5]. From M 2 and Y we reconstruct the barred momenta, as for the Born 
kinematics case. Then we reconstruct immediately 

k e = —— , k e = k e , (2-15) 

and then 

fc 3 = vk e + (1 - x - v)k @ + kr, ( 2 -16) 

where kx has only transverse components. Its magnitude is determined by the on shell 
condition fc| = 0, which yields 

k\ = 2k B ■ k e (l - x- v)v (2.17) 
and its azimuth is 4>. We then construct the vectors 

K = k B + k e -k 3 , K = xk (B + k e , (2.18) 

and the inverse boost 

\~ 1 (TC K\ — a 1{K + K)»{K + K) U 2K,K V 

%u {K, K) - (K + Kf K 2 ' ( ' 

from which we can compute the leptons' momenta 

k r = A- 1 (K,K)k r , r = l,2. (2.20) 

The real-emission phase space can be expressed in a factorized form in terms of the under- 
lying Born kinematics phase space and of the radiation variables 

d* 3 = d*2d$rad, (2-21) 

with 

^rad = ^ ^ dv ^ 9{v) 9 (l - 9(x(l -x))0(x-x e ) (2.22) 
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and 

= - dM 2 dY d cos 61. (2.23) 

6 167T 

The kinematic variables corresponding to the e collinear direction are reconstructed in full 
analogy. Observe that the underlying-Born variables and the radiation variables depend in 
general upon the collinear region that we are considering. In the present case, while M, x 
and <j) are obviously independent of the region we are considering, Y, 6\ and v do depend 
upon it. In order to avoid a too heavy notation, we have refrained from appending e or 
e indices to the underlying Born and radiation variables. When necessary, we will put a 
]@ "context" bracket around a formula, meaning that the underlying Born and radiation 
variables inside it should refer to the © direction. 

2.2 Cross sections 

We have used the helicity amplitude method of refs. [11, 12] in order to compute the cross 
sections including the vector-boson decay products. For the VF-boson propagator we have 
taken 

q 2 -M 2 r + iT w M w 
and for the Z/7-boson propagators, multiplied by the corresponding couplings, 



(2.24) 



-9»v + W M I , ziEL (2 25) 

where gi, g q are the lepton and quark couplings to the Z (for given helicities), and ej, e q 
are their electric charges. 

Following ref. [5], we introduce the Born B qq and the real-emission cross sections ~R- qq , g , 
Ttgq, q and TZqg t q, that represent the contributions for quark-antiquark, gluon-antiquark 
and quark-gluon initiating processes. Notice that the flavour of the outgoing particle in 
the subscript of 1Z is also taken to be incoming. In the case of Z production, q and q 
are conjugate in flavour. For production, because of flavour mixing, q and q may 
refer to different flavour species. We thus assume that, in general, q and q may both 
represent any flavour, but, in general, if q is a quark, q is an antiquark, and viceversa. 
B and 1Z are obtained by taking the absolute value squared of the corresponding helicity 
amplitude, summing over the helicities and colours of the outgoing particles, averaging over 
the helicities and colour of the initial partons, and multiplying by the flux factor l/(2s) 
(see eq. (2.27)). The soft-virtual term in the CS approach is given by (see eq. (2.107) in 
ref. [5]) 

Vqq = ^C F B qq . (2.26) 

Defining 

s = (fe ffi + k e ) 2 , u = (k & -k 3 ) 2 = -sv, t= {k e - k 3 ) 2 = -(l-x-v)s, (2.27) 
the CS subtraction terms are given by 



11,9 



■ h g 2 C F I ^— - (1 + x) I B qq (M, Y, SO 



(2.28) 
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for gluon radiation from a qq initial-state, and 

1 



--39,9 



-2 gj T F {1 - 2 x (1 - x)} B M -(M , F, 5,) 



(2.29) 



for the gq. Analogous formulae apply for the qq and the qg counterterms in the e collinear 
direction. 

The collinear remnants are given by 



As 



2 (1-z) 
log 



1 - z 



(I - z) 2 
(l + z)log^ tL + (l- z ) 



+( ^_ 5 W_, )+ (I±£! 



log- 



:i-z) 



1 M 

+ log — 



M 2 ' 
log— T 

21 



[B M -(M,y,e,)] ffl ,(2.30) 



+ 2z(l-z)^[B M -(M,y,0,)] e - 

(2.31) 



The *2,e notation, according to ref. [5], represents the set of variables 

*2,e = {x<B,x e ,z,ki,k 2 }, zx m K s> + x e K e = k l + k 2 . (2.32) 

We also associate an underlying Born configuration $> 2 to the $2,® kinematics, defined by 

k B = zx^K^, k e = x e K e , h = h, k 2 = k 2 . (2.33) 

The other two collinear remnants, Ge l ' 9 (&2,e) and ^| 9 ' <? (^ , 2,e)) are equal to ^® 9 ' 9 (^2,e) 
and $f' 9 (#2, e ) respectively, with [B, ff (M, Y, 0j)] replaced by [B qq {M, Y, 0j)] . We then 
introduce the notation 1?, V", i?, C, G, to stand for B, V, 1Z, C, Q, each multiplied by its 
appropriate parton densities. The differential cross section, multiplied by some infrared 
safe observable O, can then be written as 

+ J dd> 3 {^,- fl (* 3 )0(*3) " C®- fl (* 3 ) [O(* 2 )] - C^ ff (* 3 ) [0(* 2 )] e } 
+ / <**3 {^,-,(* 3 )0(*3) " <W* 3 ) [O(*2)] } 
+ | d* 3 {W*3)0(*3)-C, ff) ,-(*3) [0(*2)] Q } 

+ J d* 2 , e [Gf • 9 (* 2>e ) + Gf' «(*2, e )] 0(*2, e )} ■ (2-34) 
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3. POWHEG implementation 



The starting point of a POWHEG implementation is the inclusive cross section at fixed 
underlying-Born flavour and kinematics. For the soft-virtual and Born contributions the 
underlying Born kinematics is obviously given by the Born kinematics itself. For the 
collinear remnant, for example, in the © direction (see eq. 2.32) the underlying Born kine- 
matics is given by 

*2 = {zX(B,x e ,k 1 ,k 2 } ■ (3.1) 

For the CS counterterms, the underlying Born kinematics is given by the corresponding 
$2 variables defined in eqs. (2.11) and (2.12). In order to assign an underlying Born 
kinematics to the real term, one has to decompose it into contributions that are singular 
in only one kinematic region. Since R gq , q and R qg>q are only singular in the e and e 
direction respectively, we assign their underlying Born to be the same of the corresponding 
CS subtraction term. For R qq , g , on the other hand, we separate: 

C® 

Rqq,g = Rqq,g + Rqq,g-> ^fq,g = Rm,9~F!W_ i'ne_ ' $.2) 

^qq,g ^ ^qq,g 

and assign to R® q ^ g the same underlying Born kinematics of the corresponding CS coun- 
terterm C®q g . The underlying Born flavour, on the other hand, is always qq in the notation 
we have adopted. 



3.1 Generation of the Born variables 

The primary ingredient for a POWHEG implementation is the B function, that is the inclusive 
cross section at fixed underlying Born variables. In our case, it is given by 

B = ^B qq , (3.3) 
qq 

B qq = B qq (* 2 ) + V qq {* 2 ) + W [dt> iad {Rf q , g (* 3 ) ~ C®- 9 (* 3 )}] J 2= * 2 

© J 

+ J {iW*3) - C OT - i(? (*3)}]J 2= * 2 +J [^rad {^W^) " <W* 3 )}] 

+ f d jL [Gr(* 2 , e ) +Gr($ 2 , ffi )] #2 =* 2 [Gf ^(* 2 , e ) +G|^(* 2 , e )]*^ * 2 

Jxm z Jxq z 



#2=*2 

e 



(3.4) 

The radiation variables 'I'rad are parametrized in terms of three variables that span the unit 
cube, X ra d = {-^rad> -^rad' ^radi' wnue the z variable is parametrized in term of a single 



-7- 



variable X^ d that ranges between and 1. We then define the B function 



B qq = B qq {* 2 )+V qq (3> 2 ) + Y, 



± 



rad 



+ 



+ 



+ 



1 dz 



{^(*3)-C|- )9 (* 3 )} 
d$rad 



#2=*2 
J © 



+ 



a) 



{Gr(* 2 , e ) + Gf^(* 2 , e )} 



rad 

1 



(i) 

rad 



{G^(* 2 , e ) + G|^(* 2)e )} 



9-^rad 
#2 = ^2 

#2 = *2 



{iW*s) " <W* 3 )} 



*2 = *2 



(3.5) 



so that defining B = ^ ^ we have 



^ -^rad B . 



(3.6) 



In practice, the B function is integrated numerically over all <& 2 , -^rad integration variables, 
using an integration program that can generate the set of kinematic variables 3> 2 , A^ r ad, 
with a probability proportional to <i<fr 2 <i 3 X ra d B in the d& 2 d 3 X rax i kinematic cell (see, for 
example, refs. [13, 14]). Once the <& 2 ,X ra d point is generated, the flavour qq is chosen with 
a probability proportional to the value of B qq at that specific $ 2 , X ra d point. At this stage, 
the radiation variables are disregarded, and only the underlying Born ones are kept. This 
corresponds to integrate over the radiation variables. 



3.2 Generation of the radiation variables 

Radiation kinematics is instead generated using the POWHEG Sudakov form factor 

a^(* 2 , Pt )=;qai, 



(3.7) 



@ 



where 

Af(* 2 ,p T ) = ex P |- 

Af(* 2 ,p T ) = exp j- 



:/ 
:/ 



rad ' 



fl& )fl (* 3 ) + iW*3) 



S M -(*2) 



rad ' 



R^ gjQs) + R qg , q (*3) 



(fc T (* 3 ) - Pt ) 



9 (fcr(* 3 ) - Pt) 



#2=*2 



#2=*2 



(3.8) 
(3.9) 



The function &t(* 3 ) measures the hardness of radiation in the real event. It is required 
to be of the order of the transverse momentum of the radiation in the collinear limit, and 
to become equal to it in the soft-collinear limit. In principle, the choice of fcr(^ 3 ) can 
differ in the two singular regions (e and e) that we are considering. The choice adopted 
in the Examples section of ref. [5] had in fact this feature. We have found, however, that 
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for practical reasons 3 it is better to adopt a different choice, namely to take kri&s) to 
coincide with that of eqs. (2.16) and (2.17). 

The generation of radiation is performed individually for A| 9 and A^, and the highest 
generated kx is retained. The upper bounding function for the application of the veto 
method is chosen to be 4 

B q q ~ M2 ^ 2V 1-X-V' { ' 

and the analogous one for the e direction. The procedure used to generate radiation events 
according to this upper bounding function is described in Appendix A. 

3.3 Born zeros 

In case the Born cross section vanishes in particular kinematics points, a problem arises 
in the POWHEG expression for the Sudakov form factor (3.8) and (3.9). It happens, in 
fact, that although B vanishes, B may differ from zero. Born kinematics configurations 
with a vanishing Born cross section may thus be generated and, at the stage of radiation 
generation, one would find very large ratios of the real-emission cross section over the 
Born cross section. It would thus prove difficult to find a reasonable upper bound for this 
ratio. If one tries to neglect the problem, radiation events with a vanishing underlying 
Born configuration would never be generated. We observe that, in the limit of small 
hardness parameter, the real cross also exhibit the same vanishing behaviour of the Born 
cross section. Loosely speaking, the problem arises when the distance of the underlying 
Born configuration from the zero configuration is smaller than the distance of the real 
emission cross section from the singular (i.e. zero hardness) configuration. In order to solve 
this problem, in a completely general way, we further decompose the real cross section 
contribution as (we use the notation of ref. [5]) 

R a r = R a r ,s + R a r ,r^ (3 n ) 

where 

R a r ,s = R a r R a r ,r = R a r (3^) 

Z + H' Z + H v ; 

The suffixes s and r stand for "singular" and "regular" respectively, and Z is a function 
of the kinematics that vanishes like the Born cross section, evaluated at the underlying 
Born kinematics of the given term. H is the hardness of radiation and it must vanish for 
vanishing transverse momentum of the radiation. The simplest possible choice would be 

k 2 

Z = B-^^, H = k T , (3.13) 



The choice discussed in [5] is k T — M (1 — x)v, and is such that k T is always bound to be smaller than 
M 2 . Since the factorization and renormalization scales are taken equal to kr, for vector-boson production 
at transverse momenta much larger than the vector-boson mass the coupling does not properly decrease. 

4 This upper bounding function differs from the ones of eqs. (7.163)-(7.166) in ref. [5], but is in fact 
equivalent to the bound of cq. (7.234) in the same reference, once the change of variables £ = 1 — x, 
y = (l — 2v — x)/(l — x) is performed, and the different definitions of d$ ra d are properly taken into account. 
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where kx is some definition of the transverse momentum of the radiation. Notice now that 
R ar ' s vanishes as fast as the Born term when its underlying Born kinematics approaches 
the Born zero. It can thus be used in the expression for the Sudakov form factor (eqs. (3.8) 
and (3.9)) without problems. The R ar > r is instead non-vanishing, but, on the other hand, 
it does not have collinear or soft singularities because of the H factor, and thus it can be 
computed directly, without any Sudakov form factor. In the case of W production, the 
Born zero is associated to Oi = if q is an antiquark, and Oi = ir if it is a quark. We choose 
then 

Z = M 2 (l + s q cos Oif , H = k\, (3.14) 

with k\ given by formula (2.17) and the factor s q equals 1 for quark, and —1 for antiquark. 
The angle Q\ is chosen according to the ® parametrization (for R®) or the e parametrization 
(for R G ) of the real-emission phase space. 

In addition, all the R ar terms in eq. (3.3) are replaced by the corresponding R ar ' s and 
the R ar ' r terms are generated in a way similar to what was done for eq. (3.5). In other 
words one defines 



d<S> 



rad 



+ 



dX rad 



R^ r r*,i 



dX, 



rad 



R r - 



#2 = ^2 



3>2=*2 



rad 



+ 



5$rad 



dX, 



rad 



^qq^y^Zj 



"19,1 



*2=*2 



<I>2 = *2 



, (3.15) 



and integrates over the whole <&2,^rad phase space with the same method used for B. In 
order to generate an event, one chooses B or B r , with a probability proportional to their 
respective total integral. In case B r is chosen, one generates a kinematic configuration 
according to it. This kinematic configuration is a full 3-body configuration. The flavour qq 
is chosen with a probability proportional the the value of B r g for the particular kinematic 
point that has been generated, and the event is sent to the output. In case B is chosen, a 
kinematic configuration and an underlying Born flavour is chosen in the same way. 



4. Results 

The MCONLO program provides an implementation of vector-boson production at the NLO 
level in a shower Monte Carlo framework. It should therefore be comparable to our calcula- 
tion, and we thus begin by comparing MC@NL0 and POWHEG distributions. In this comparison, 
the POWHEG code is interfaced to HERWIG [15, 16], in order to minimize differences due to 
the subsequent shower in the two approaches. We choose as our default parton-density 
functions the CTEQ6M [17] package, and the corresponding value of A QCD . The factor- 
ization and renormalization scales are taken equal to My + (p^f) 2 in the calculation of the 
B function, where V = W or Z. In the generation of radiation, the factorization and 
renormalization scales are taken equal to the transverse momentum of the vector boson V . 
We also account properly for the heavy-flavour thresholds, when the transverse momentum 
of the vector boson approaches the bottom and charm quark threshold. That is to say, 
when the renormalization scale crosses a heavy-flavour mass threshold, the QCD evolution 
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of the running coupling is accordingly changed to the new number of active flavours. The 
other relevant parameters for our calculation are 



M z (GeV) 


T z (GeV) 


M w (GeV) 


rv (GeV) 


sin 2 9§ 




91.188 


2.49 


80.419 


2.124 


0.23113 


127.934 



The above values of masses and widths are used in eqs. (2.24) and (2.25). The W and Z 
couplings are given by 

3 = — 4rfF ' 9l/ q = . a J ^ff T^ l/q) - q lf sin 2 6<§ 

sin Vyy 1 sin 6yy cos 9^} L ' H 

where I / q denotes the given left or right component of a lepton or a quark. For W produc- 
tion we used the following absolute values for the CKM matrix elements 



ud 


us 


ub 


cd 


cs 


do 


td 


ts 


tb 


0.9748 


0.2225 


0.0036 


0.2225 


0.9740 


0.041 


0.009 


0.0405 


0.9992 



In all figures shown in the following we do not impose any acceptance cut. 
4.1 Z production at the Tevatron 

In fig. 1 we show a comparison of the lepton transverse momentum and rapidity, and 
of the transverse momentum of the reconstructed lepton-antilepton pair at the Tevatron. 
We notice a larger cross section in POWHEG, when the Z transverse momentum becomes 
large. This is not unexpected, since for large momenta the POWHEG result is larger than 
the standard NLO result by a factor B/B (this feature has also some impact upon the 
transverse- momentum distribution of the lepton). Once this fact is accounted for, the 
transverse- momentum distribution of the Z is in fair agreement, although we find observ- 
able shape differences at low transverse momenta. We also notice a peak at pr = in the 
MC@NL0 distribution, that is not present in the POWHEG result. We expect this distribution 
to be affected by low transverse-momentum power-suppressed effects. In fact, the peak at 
zero transverse momentum in MC@NL0 disappears if the primordial transverse momentum of 
the partons (the PTRMS variable in HERWIG) is set to a non-zero value. In fig. 2 we compare 
the rapidity distribution of the reconstructed Z, its invariant mass, the azimuthal distance 
of the e + e~ pair coming from Z decays, and the transverse momentum of the radiated jet 
at the Tevatron. The jet is defined using the SISCONE algorithm [18] as implemented in 
the FAST JET package [19], using R = 0.7. We find again fair agreement. 

In ref . [20] , a discrepancy was found in the rapidity distribution of the hardest radiated 
jet as computed in MC@NL0 and ALPGEN, for the case of top pair production at the Tevatron. 
The MC@NL0 calculation shows there a dip at zero rapidity, not present in ALPGEN. In fact, 
the POWHEG calculation of this quantity does not display any dip. We thus examine the 
transverse momentum of the radiated jet in this case. Furthermore, we also plot the rapidity 
difference between the Z and the hardest radiated jet. The results are displayed in fig. 3. 
We have chosen different cuts for the minimum transverse momentum of the radiated jet, 
i.e. 10, 20, 40, 60 and 80 GeV. We observe noticeable differences in the rapidity distribution 
of the hardest jet in the two approaches. The MC@NL0 result displays a dip at zero yj e t — Vz- 



e= v / 47ra cm (M z ), (4.1) 
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Figure 1: Comparison between POWHEG and MCONLO results for the transverse momentum and 
rapidity of the lepton coming from the decay of the Z boson, and for the transverse momentum 
of the Z, as reconstructed from its decay products. The lepton-rapidity asymmetry is also shown. 
Plots done for the Tevatron pp collider. 



4.2 Z production at the LHC 

Similar results are reported for the LHC in fig. 4 through 6. We notice less pronounced 
differences (with respect to the Tevatron case) in the pt spectrum of the Z boson. The 
discrepancy in the j/j et distribution is still evident, although the dip is barely noticeable in 
this case. 

The same set of plots are also shown for a PYTHIA-P0WHEG comparison in fig. 7 through 10. 
In this case the POWHEG code was interfaced with PYTHIA. Photon radiation from final-state 
leptons was switched off (MSTJ(41)=3), in order to simplify the analysis. Furthermore, 
the new transverse- momentum ordered shower was used (i.e. the PYEVNW routine), since 
transverse-momentum ordering should be more appropriate in conjunction with POWHEG. 
In the plots, the PYTHIA output is normalized to the POWHEG total cross section. From 
fig. 7 through 8, we can see a remarkable agreement between the two calculations for the 
Tevatron results, the only visible discrepancy being given by the transverse-momentum 
distribution of the Z boson at small transverse momenta. We also notice that, unlike the 
case of the MC0NL0-P0WHEG comparison, the transverse-momentum distribution of the Z is 
slightly harder in PYTHIA than in POWHEG. The rapidity distributions of the hardest jet are 
also in remarkable agreement. 
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Figure 2: Comparison between POWHEG and MC@NL0 for the reconstructed Z rapidity, its invariant 
mass, the lepton-pair azimuthal distance and the transverse momentum of the reconstructed jet, 
above a 10 GeV minimum value. 




In fig. 10 through 11, we carry out the same comparison in the LHC case. We notice 
here few important differences in the rapidity distribution of the Z boson, and, probably 
related to that, of the electron, the PYTHIA distribution being natter in the central region. 
Both MC@NL0 and POWHEG do not show this feature. As already pointed out in ref. [5], the 
generation of vector bosons in PYTHIA is not very different from the POWHEG generation. 
Radiation is generated with a very similar method [21, 22]. There are however differences. 
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In PYTHIA the Born inclusive cross section is used rather than our B function. Further- 
more, our choice of scales is constrained by the requirement of next-to-leading logarithmic 
accuracy in the Sudakov form factor. The discrepancy in the transverse-momentum dis- 
tribution of the Z may be due to different requirements for the choice of the scale in the 
generation of radiation in the two algorithms. The discrepancy in the rapidity distribution 
may be due to the lack of NLO corrections in PYTHIA, i.e. to the use of the Born cross 
section (rather than the B function) and LO parton densities. In fact, in fig. 3 of ref. [23], 
a comparison in the rapidity distribution of the Z at LO, NLO and NNLO, is shown for the 
LHC. One can notice from that figure that there is a difference in the LO and NLO shape 
of the distribution, the former being flatter. In order to elucidate this point, we show in 
fig. 13 the rapidity distribution of the Z boson computed at fixed order in QCD, at LO 
and NLO. With the LO calculation, we also show the result obtained using the same LO 
parton-distribution function (pdf) set used in PYTHIA, that is CTEQ5L. The figure leads 
to the conclusion that the use of the LO parton-density set CTEQ5L is the primary cause 
of this shape difference. We find, in fact, no difference in shape between the LO and NLO 
result if the same pdf set is used instead. We thus conclude that also the effect observed in 
fig. 3 of ref. [23] is due to the use of a LO parton-density set together with the LO result. 5 
The predictions for the transverse-momentum distribution of the Z boson are sum- 



Some authors do prefer to use LO parton-density functions in LO calculations, although, in our opinion, 
there are no compelling reasons to do so. 




Figure 5: Same as fig. 2 for the LHC at 14 TeV. 




Figure 6: Same as fig. 3 at the LHC at 14 TeV. 



marized in fig. 14, in comparison with data from ref. [24], at = 1960 GeV and from 
refs. [25, 26, 27] at v 7 ^ = 1800 GeV. The POWHEG+HERWIG and the MC0NL0 output are 
obtained with an intrinsic transverse momentum of the incoming partons equal to 2.5 GeV 
(HERWIG's PTRMS parameter). Both data and predictions are normalized to 1. The dif- 
ference in the shape of the distributions at 1960 and 1800 GeV are only minimal. We 
see that P0WHEG with PYTHIA is in remarkable agreement with the MC0ML0 result. On the 
other hand, standalone PYTHIA is closer to the output of P0WHEG with HERWIG. In all cases, 
the agreement with data is not optimal. It is thus clear that this distributions is sensi- 
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tive to long distance effects like hadronization and transverse-momentum smearing, and 
good agreement with data may only achieved by suitable tuning of the non-perturbative 
parameters of the shower Monte Carlo. 




1 2 3 50 100 150 200 250 300 



A*(e + e") pf 
Figure 8: Same as fig. 2 for a PYTHIA and PDWHEG comparison at the Tevatron. 




yjet yjet Yz 



Figure 9: Same as fig. 3 for a PYTHIA and PDWHEG comparison at the Tevatron. 
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Figure 12: Same as fig. 3 for a PYTHIA and PDWHEG comparison at the LHC. 
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Figure 13: Rapidity distribution for the Z boson, computed at fixed order at LO and NLO. For 
the LO result, both the CTEQ6M and the CTEQ5L parton-density set were used. The plots are 
normalized to the NLO total cross section. 
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4.3 Hardest-jet rapidity distribution 

The discrepancy of POWHEG and MC@NL0 in the rapidity distribution of the hardest jet 
deserves further discussion. 6 In ref. [20], only the rapidity distribution of the hardest jet in 
tt production was considered, and a dip was found there, in the case of top-pair production 
at Tevatron energies. In the present case we found no dip in the rapidity distribution of 
the hardest jet in V production (see fig. 3). We found instead a dip in the distribution in 
the rapidity difference between the jet and the vector boson. It is reasonable to assume 
that a dip in the rapidity distribution of the jet may be inherited from the dip in the 
rapidity difference, if the kinematics production regime is forced to be central, like in the 
case of top-pair production at the Tevatron. We thus also reconsider Z pair production 
and production at the Tevatron, and compare POWHEG and MC@NL0 results for the rapidity 
distribution of the hardest jet, and for the distribution in the rapidity difference. The results 
are shown in figs. 15 and 16. From fig. 15 we see that the dip present in the j/j e t — Vtt 
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Figure 15: Rapidity distribution of the hardest jet and of the rapidity difference between the 
hardest jet and the tt system at Tevatron energies. 



distribution is even deeper than the dip observed in the yj e t distribution. Furthermore, in 
fig. 16, we see no particular features in the yjet distribution. The y- ]e t — yzz distribution 
displays instead a tiny tower and a dip, depending upon the transverse-momentum cut on 
the jet. A deeper study of these features was performed in ref. [20], for tt production. It 
was shown there that the HERWIG Monte Carlo displays an even stronger dip than MC@NL0. 
The MC@NL0 generator provides more events that partially fill the dip, thus correcting the 
NLO inaccuracies of the shower Monte Carlo. It is presumably a NNLO (next-to-next-to 
leading) mismatch between the twos that generates these features. On the other hand, the 
POWHEG program, as well as matrix-element generators, generate themselves the full NLO 
result, and thus are not sensitive to this feature of HERWIG. We also stress that these features 
do not mean that HERWIG is inaccurate at the LO level, or that MC@NL0 is inaccurate at 
the NLO. A shower Monte Carlo is accurate in the radiation of the hardest jet only in the 



6 The distribution in the pseudorapidity difference of the hardest jet with respect to the vector boson 
was considered in ref. [28], in the context of a comparison of several matrix-element programs. Although 
noticeable differences are found among the generators considered there, none of them exhibit a dip at zero 
pseudorapidity. 
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Figure 16: Rapidity distribution of the hardest jet and of the rapidity difference between the 
hardest jet and the ZZ system at Tevatron energies. 

collinear regions. Furthermore, the dip in the MC@NL0 result is compatible with an effect 
beyond NLO. 

4.4 W production at the Tevatron and LHC 

All results presented so far are relative to Z boson production. In the case of W production 
we find similar features and the comparison between MCONLO and PYTHIA presents very sim- 
ilar characteristics. For the sake of completeness, we present in fig. 17 through 34 plots of 
observables for W~ production at the Tevatron, and W~ and W + production at the LHC, 
comparing again the POWHEG output with MCSNLO and PYTHIA, and the observables for W + 
production at the LHC. We find again that MCOMLO displays dips in the rapidity distribu- 
tion of the hardest jet at Tevatron energy. The comparison of the transverse- momentum 
distribution of the W shows the same differences found in the Z case. Furthermore, the 
rapidity distribution of the at the LHC differs in PYTHIA, showing a very marked 
difference in the W + case (see fig. 33), probably (as in the Z case) a consequence of the 
different pdf set. 
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Figure 17: Comparison of POWHEG and MCONLO results for the transverse momentum and rapidity 
of the lepton coming from the decay of the W~ boson and for the transverse momentum of the 
W~ , as reconstructed from its decay product. 
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Figure 18: Comparison of PDWHEG and MC0NL0 for the reconstructed W~ rapidity, its invariant 
mass, the lepton-pair azimuthal distance and the transverse momentum of the reconstructed jet, 
above a 10 GeV minimum value. 
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Figure 19: Rapidity distribution of the hardest jet with different transverse-momentum cuts, and 
the rapidity distance between the hardest jet and the reconstructed W~ boson for POWHEG and 
MC0NL0 at the Tevatron. 
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5. Conclusions 

In this paper we have reported on a complete implementation of vector-boson production 
at NLO in the POWHEG framework. The calculation was performed within the Catani- 
Seymour [6] dipole approach, and thus this is the first POWHEG implementation within the 
Catani-Seymour framework at a hadronic collider. We have found that, at variance with 
what was proposed in sec. 7.3 of ref. [5], it is better to define the transverse momentum as 



- 30 - 




-4 -2 2 4 50 60 70 80 90 100 110 

M(i/e + ) (GeV) 




1 2 3 50 100 150 200 250 300 

Figure 33: Same as fig. 27 for W + production at the LHC. 
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Figure 34: Same as fig. 28 for W + production at the LHC. 

the true transverse momentum for the initial-state singular region. Furthermore, we have 
shown how to perform a PDWHEG implementation when the Born term vanishes. 

The results of our work have been compare extensively with MC@NL0 and PYTHIA. The 
PYTHIA result, rescaled to the full NLO cross section, is in good agreement with P0WHEG, 
except for differences in the rapidity distribution of the vector boson, that may be ascribed 
to the use of a LO parton density in PYTHIA. The MC@NL0 result is in fair agreement 
with PDWHEG, except for the distribution of the hardest jet in the process, the MC@NL0 
distribution being generally wider. Furthermore, we have also examined the distributions 
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in the difference of the hardest jet and the vector-boson rapidity. We have found that these 
distributions exhibit dips at zero rapidity. We have also examined analogous distributions 
for ZZ and ti production, and again found dips for these distributions, that seem to be a 
general feature of the MCONLO approach. We also remark that no other approaches show- 
dips of this kind [28] . 

The computer code for the PDWHEG implementations presented here is available, to- 
gether with the manual, at the site 

http : //moby .mib . inf n. it/~nason/PDWHEG . 
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A. Upper bounding function 

We call Ajj(pt) the Sudakov form factor obtained with the upper bounding function of 
eq. (3.10). Using the definitions of eqs. (2.22) and (2.17) 

d^rad = ^ ^ dv J 9{v) 9 (l - 9(x(l - x)) 9(x - x e ) (A.l) 

M 2 

kl = (l-x-v)v , (A.2) 

x 

we write 

logAuiti) [Ux f 1 '* a s {k 2 ) x 2 2 2 

-N =L V 2 L dv 2v i-x-y ^-Pr) 
f 1 dx f l - x , a s (kl) M 2 , 2 2X 

f°° dk 2 a s (k 2 ) f 1 , f 1 dx n , . o J M 2 / N , 2 \ 

where, for ease of notation, we have dropped the © and qq labels on N and x. We perform 
the x integration using the 5 function 

Notice that x < 1, and 

/ M 2 ^(l-t;)\ / 4 + M 2 \ 1 , A ^ 

The only remaining condition on x is x ^ x. We thus get 

logAf/^) f°° dk 2 a s (kl) f 1 dv (M 2 v{l-v) 



r dk 2 T a s (k 2 ) f l dv J M 2 v{l-v) _\ 
L k 2 2 J k 2 T /M 2 + V \k 2 +M 2 v X J- ( Dj 
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We must find the conditions implied by the theta function upon v. For 

M 2 (l-x) 2 



k\ < fe|i max — 



Ax 



(A.6) 



the function is satisfied if V- < v < v+, where 



l-x± 



- x) 



2 -4: 



M 2 



(A.7) 



We thus have 



log Ag(p|) 
-N 



y*2 



2 



log 



M 2 
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(A.8) 



The k\ integral is still too complex to be performed analytically. We thus resort another 
time to the veto method, by finding an upper bound to the integrand. We have 



k 2 

-% + V 



M 2 



kl k 2+ ^k 2 T + k 2 
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M 7 



M 2 (l + x) 
Axk\ 



(A.9) 



We thus define 



M 2 (l+aQ 2 
Axk\ ' 



and introduce a new Sudakov form factor 



(A.10) 



logA^O 2 .) _ f^Tmax dk 2 T a v (k 2 T ) q 
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Pt 



l0g P"' 

IX rp 



(A.ll) 



where au(k^) has the form of the one-loop running coupling constant 



k 2 ' 
&l0g-rf 



(A.12) 



and is required to satisfy the bound au(k^) ^ d s {k 2 -) in the allowed range for k\. The 
integral in eq. (A.ll) is now easily performed, and we get 



&u(Pt) = ex P * 



N 
2b 



log -2- lot 
1V U 



log^f- ^ k 
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log |f Pt 



(A.13) 



The generation of the radiation variables is then performed starting with Au(p 2 ^), using 
the veto procedure to obtain the Au(p 2 ^) distribution. Further vetoing is then used to 
obtain the correct R/B generated distribution. 
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